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Abstract 

In this work we deal with the optimal design and optimal control of structures under- 
going large rotations. In other words, we show how to find the corresponding initial 
configuration and the corresponding set of multiple load parameters in order to re- 
cover a desired deformed configuration or some desirable features of the deformed 
configuration as specified more precisely by the objective or cost function. The model 
problem chosen to illustrate the proposed optimal design and optimal control method- 
ologies is the one of geometrically exact beam. First, we present a non-standard for- 
mulation of the optimal design and optimal control problems, relying on the method 
of Lagrange multipliers in order to make the mechanics state variables independent 
from either design or control variables and thus provide the most general basis for de- 
veloping the best possible solution procedure. Two different solution procedures are 
then explored, one based on the diffuse approximation of response function and gra- 
dient method and the other one based on genetic algorithm. A number of numerical 
examples are given in order to illustrate both the advantages and potential drawbacks 
of each of the presented procedures. 
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1 Introduction 

Modern structures should often be designed to withstand very large displacements and 
rotations and remain fully operational. Moreover, the construction phase has also to 
be mastered and placed under control, with trying to precisely guide the large motion 
of a particular component of the total structural assembly. Ever increasing demands 
to achieve more economical design and construction thus require that the problems of 
this kind be placed on a sound theoretical and computational basis, such as the one 
explored in this work. Namely, optimization methods can be called upon to guide 
the design procedure and achieve desired reduction of mechanical and/or geometric 
properties. Similarly, the control methods are employed to provide an estimate of the 
loads and the minimal effort in placing the structure, or its component, directly into 
an optimal (desired) shape. Either of these tasks, optimal design or optimal control, 
can formally be presented as the minimization of the chosen cost or objective function 
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specifying precisely the desired goal. The main difference between two procedures 
concerns the choice of the variables defining the cost function: the design variables 
are typically related to the mechanical properties (e.g. Young's modulus) or geometry 
of the structure (e.g. particular coordinates in the initial configuration), whereas the 
control variables are related to the actions (e.g. forces) applied on the structure in order 
to place it into desired position. Rather then insisting upon this difference and treating 
the optimal design and optimal control in quite different manners (as done in a number 
of traditional expositions on the subject), we focus in this work on their common 
features which allow a unified presentation of these two problems and the development 
of a novel solution procedure applicable to both problems. The latter implies that the 
nonlinear mechanics model under consideration of geometrically exact beam has to 
be placed on the central stage and one should show how to fully master the variation 
in chosen system properties or loads in order to achieve the optimal goal. The main 
contributions put forward in presenting this unconventional approach can be stated as 
follows: 

i) We first present the theoretical framework for treating nonlinear structural me- 
chanics and optimization and/or control as a coupled (nonlinear) problem; In 
any such problem of optimization or control, the nonlinear mechanics equilib- 
rium equations are reduced to a mere constraint with respect to the admissibil- 
ity of a given state of the structure, i.e. its displacements and rotations. By 
using the classical method of Lagrange multipliers (e.g. see [fT9l or ||29l the 
mechanics equilibrium equations can be promoted from constraint to a one the 
governing equations to be solved in a coupled problem of this kind, and the in- 
trinsic dependence on state variables (displacements and rotations) with respect 
to optimal design or control variables can be eliminated turning all variables 
into the independent variables. For clarity, this idea is also developed within 
the framework of a discrete, finite-element-based approximation, thus provid- 
ing the finite element model including the degrees of freedom pertinent not only 
to displacements and rotations, but also to optimal design and/or control vari- 
ables. A detailed development is presented for the chosen model problem of 3D 
geometrically exact beam (e.g. see 11271 or ifTOlO . We note in passing that the 
proposed approach is quite opposite from the traditional ones (e.g. see ifTTl for 
a very recent review), where the two fields directly concerned by the coupled 
problems of this kind, nonlinear mechanics on one side and optimal design or 
control on the other, are studied and developed separately and resulting solution 
procedures for one or another are applied in a sequential manner. One typically 
employs two different computer programs, one for mechanics and another for 
optimization or control (e.g. the toolbox of MATLAB), so that the communi- 
cation requirements are reduced to a bear minimum: so-called sensitivity (e.g. 
11301 or I1251D for optimization code, and design or control variables for the fi- 
nite element code for mechanics. It is clear that such a traditional approach to 
analysis and design or control will largely sacrifice the computational efficiency 
for the cases of practical interest where both the cost function and mechanics 
problem are nonlinear and require iterative procedures to be solved. 

ii) The second aspects which is elaborated upon in this work pertains to an al- 
ternative method for solving a coupled problem for analysis and either design 
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or control, where two sub-problems, once brought up to the same level by the 
Lagrange multiplier method, are solved simultaneously. In other words, once 
the interdependence of the state variables, i.e. the displacements and rotations 
on one side and the design or control variables or another, is no longer en- 
forced, one can iterate simultaneously to solve for all of them. In particular, 
the sensitivity analysis needs no longer be performed separately, but it is natu- 
rally integrated as a part of the simultaneous iterative procedure. It is important 
to note that the iterative intermediate values are no longer consistent with the 
equilibrium equations constraint, except at the convergence, where basically the 
same solution is obtained as for the standard sequential solution procedure but 
with a (significantly) reduced number of iterations. In other words the simul- 
taneous and sequential solution procedures will always yield the same result 
providing the solution is unique. The problems of this kind concern the opti- 
mization and/or control of geometrically nonlinear response of structures where 
no bifurcation phenomena occur. Solving the nonlinear problems in structural 
mechanics can be quite a demanding task for the structures where the stiffness 
may differ significantly in different deformation modes (i.e. beam bending ver- 
sus stretching). Adding the control or optimization problem on top makes such 
a nonlinear problem so much more challenging. Different modifications of the 
basic solution procedure are developed and tested, including genetic algorithms 
to explore the initial phase of the solution procedure, gradient based ac- 
celeration near solution points and response surface for a part of the solution 
constructed by diffuse approximation (e.g. see 0). 

The outline of the paper is as follows. In the next section we briefly present the 
chosen model of geometrically exact 3D beam, capable of representing very large 
displacements and rotations. The theoretical formulation of the optimal control and 
optimal design of the chosen mechanical model is presented in Section [3j along with 
their discrete approximations constructed by the finite element method given in Sec- 
tion HI The proposed solution procedure is described in detail in Section [5l Several 
numerical examples are presented in Section[6]in order to illustrate the performance of 
various algorithms used in computations. Some closing remarks are stated in Section 

m 

2 Model problem: geometrically exact 3D beam 

In this section, we briefly review the governing equations of the chosen model problem 
of a structure undergoing large rotations, a 3D curved beam. For a more thorough 
discussion of the chosen model we refer to E71 . Il26l or |[T3~1 . among others. We 
assume that the initial configuration of the beam is internal force free and that it can 
be described by a 3D position vector <p (s) identifying the position of each point of the 
neutral fiber (an inextensible fiber in pure bending) and the corresponding placement 
of the cross-section of the rod, which is carried out by choosing a local ortho-normal 
triad of vectors. The vector triad of this kind can be obtained by simply rotating the 
global triad by an orthogonal tensor.For a usual choice of so-called normal coordinates 
with the first vector of the triad being orthogonal to the cross-section and the remaining 



3 




0o = Oo> Ao)\ 
£ 1 



\ /0-(V:A) 
] 



Figure 1 : Initial and deformed configuration of the 3D geometrically exact beam. 

two placed in the plane of the cross section, this orthogonal tensor becomes a known 
function of the initial configuration, A (s). For the case of a curved beam studied 
here, V is chosen as the arc length parameter. 

By applying the external loading f(i), parameterized by pseudo-time 't' (where 
'pseudo' implies that the inertia effects are neglected) we obtain the beam deformed 
configuration defined by the position vector cp(s, t) and the orthogonal tensor A(s, t). 
The latter is an accordance with the usual kinematic hypothesis that the cross-section 
of the beam would not deform, which, along with the hypothesis that the first vector in 
the triad remains orthogonal to it (with other two within the plane of the cross-section) 
fully determines A(s, t). See Figure 1 where the initial and deformed configuration of 
the beam are presented. 

In short, one can state that the configuration space of the described model of 3D 
beam can be written as 

C:={(f> = (cp,A)\<peR 3 ,Ae SO(3)} (1) 

where R 3 and 5*0(3) are spaces of 3D vectors and special orthogonal tensors, respec- 
tively. 

The main difficulty in numerical solution of the structural mechanics problems 
featuring the beams of this kind stems from the presence of SO (3) group in its config- 
uration space (e.g. see flU, tfT3l . |fTT]| for a more thorough discussion of these issues). 
In short, in performing the standard task of computing the virtual work principle or the 
consistent linearization, where a small rotation described by a skew-symmetric tensor 
5® E so(3), ought to be superposed on a large rotation described by an orthogonal 
tensor A 6 SO (3) one must first make use of exponential mapping 

A £ = Aexp\e5®}; exp[5®] = cos 50 1 + ^^-5® + 1 ~ °° S — SO ® 56 (2) 

50 50 

where 56 is the axial vector of the skew- symmetric tensor 5® i.e. 5® v = 56 x v, 

Vv G M 3 . 

The complexity of the last expression is in sharp contrast, with respect to a sim- 
ple additive update of virtual displacement field 5<p E M 3 when superposed on the 
deformed configuration ip E M 3 

<Pe = <P + £5V (3) 
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The results in (2) and (3) can be presented in an equivalent form by stating that the 
tangent space of the chosen beam model is defined by 



TC := {5(f) := (6<p, 56) \5<p G R 3 , 50 G M 3 } (4) 



The strain measures employed in this beam theory (e.g. IfTOlO can be written in 
direct tensor notation form as 

e = (5) 

for the axial and shear strains, and 

Q = A T A' ; ftv = u3 x v ; Vv G K 3 (6) 

for bending and torsional strains. In (5), (6) and subsequent equations we denote with 
superposed prime the derivative with respect to arc-length coordinate in the initial 
configuration, i.e. 

We consider the simplest case of linear elastic material model for the beam which 
allows us to express the constitutive equations in terms of stress resultants as 

n=C(e-e ); C = diag(EA, GA, GA) (8) 

m = D(w-w ); D = diag(GJ, EI, EI) (9) 

For illustration, we also consider the simplest case of a circular cross section, with 
section diameter 'ef , and 

d 2 -K _ d 4 7i T d 4 7i 

A = ; 1 = ; J= (10) 

4 ' 64 ' 32 

as the section area, moment of inertia and polar moment. 

In order to complete the description of the chosen beam model we state the equi- 
librium equations in the weak form as 

5(f)) := J(5e-n + 5u- m) ds - G ext {5<t>) = (11) 

where G ext (5(p) is the external virtual work and 5e and 5ui are the virtual strains. The 
latter can be obtained as the Gateaux derivative of the real strains in (5) and (6), by 
taking the results in (2) and (3) into consideration. In particular, this leads to 



5e{^54>) = 

= (12) 
= A T V + ex56 
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and 



= ^ A IKJ l-o d3) 

= <5e' + <50 T fi + me 

which can also be written in an equivalent form in terms of the corresponding axial 
vectors 

5u = 56' + uj x 56 (14) 
For the external load which derives from a given potential, such that 

G ext (5<j)) := Z^IW^) (15) 

one can also define the governing equilibrium equation of the chosen model problem 
of the geometrically exact elastic beam from the principle of minimum of the total 
potential energy defined according to 

11(0) := J^{( e - Co) • n + (u> - u> ) • m} ds - U ext ((j>) -> min (16) 
which implies that 



D4n(0)] = G(0;<ty) = O 
11(0) = minn(0*) <J D^D^cj))} = D^Gj^; 5<j>) > (17) 



V</>* 



In (17) above we denoted by K the second variation of the total potential energy 
or the tangent operator, which is obtained by the consistent linearization procedure 
(e.g. see 11201 ). For pertinent details of the consistent linearization we refer to [TTOl for 
quaternion parameterization of finite rotations or lfT3l for rotation vector- like parame- 
ters for finite rotations. 

An alternative case of external loading where the potential may not be defined and 
where the weak form in (1 1) is rather the starting point of the solution procedure, can 
also be encountered in applications. For example, in the case of the follower force po 
and follower moment 1 which follow the motion of a particular cross-section of the 
beam at node 'a', we can express their contribution to virtual work according to 

G ext (A a ; 5<p a , 56 a ) := 5ip a ■ A Po +56 a ■ Al (18) 

p i 

The follower force and moment will also contribute to the tangent operator accord- 
ing to 

D^G ext {K- 5<p a , 56 a ) ■ ( ^ J = 6<p a ■ PA6 a + 56 a ■ LA0 Q (19) 

where Pv = p x v; Lv = 1 x v; Vv G IR 3 . This contribution should also be taken 
into account when computing the solution to (11) and trying to ensure the quadratic 
convergence rate. 
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3 Coupled optimality problem 



The presented beam model provides an excellent basis to master the optimization 
problem as well as the control problem of geometrically nonlinear elastic structures. 
Although the former deals with geometric characteristics of the beam and the latter 
with the external loading sequence, the two problems can be formulated and solved in 
quite an equivalent manner, as shown next. 

3.1 Optimal design 

The optimal design problem addressed herein pertains to selecting the desired features 
of the mechanical model by leaving the free choice of the geometric properties of the 
beam model, e.g. the thickness variation or the chosen initial shape. This task is often 
referred to as the shape optimization. From the mathematical standpoint the shape 
optimization can be formulated as the problem of minimization of so-called objective 
or cost function «/(•), specifying the desired features. The latter is considered as a 
functional which depends not only on mechanics state variables <f> but also on design 
variables d (e.g. the beam thickness or its shape). 

The shape optimization procedure is interpreted herein as minimization of so-called 
cost or objective functional J(-), which can be written as 

J((j>(d),d) = min j(<t>(d*),d*) (20) 
G(0V);-)=o 

Contrary to the minimization of the total potential energy functional in (17), not all 
mechanical and design variables are admissible candidates, but only those for which 
the weak form of the equilibrium equations is satisfied. In other words, we now need 
to deal with a constrained minimization problem. 

The classical shape optimization procedure of solving this constrained minimiza- 
tion problem is carried out in a sequential manner, where for each iterative value of 
design variables dW, a new iterative procedure must be completed leading to 4>{d^) 
verifying the equilibrium equations. The considerable computational cost of such a 
procedure, most of it waisted on iterating to convergence on equilibrium equations 
even for non-converged values of design variables, can be drastically reduced by for- 
mulating the minimization problem in (|20l) using the method of Lagrange multipliers 
with 

max min L(f,d*; A); L(0*, d*; A) = J(</>*, d*) + G((f)*, d*, X) (21) 

VA V(0*,d*) 

In (|2T|) above A are the Lagrange multipliers inserted into the weak form of equi- 
librium equations instead of virtual displacements and rotations. In accordance to the 
results presented in the previous section we can write explicitly 

d- X) = j{X' v ■ An + X e ■ (E T n + ft T m) + A^ ■ m}ds - G ext (X) (22) 
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The main difference of (|2TI) with respect to constrained minimization problem in 
(1201) pertains to the fact that state variables <fi and design variables d are now consid- 
ered independent and they can be iterated upon (and solved for) simultaneously. 

The Kuhn-Tucker optimality condition (e.g. lfT9l ) associated with the minimization 
problem in (18) can be written as 



= [L(<f>, d; A)] = [J{<t>, d)} + [G(<f>, d; A)] (23) 
with the explicit form of the last term which can be written as 

[G(cf>, d; A)] ■ 5(f> = f t {A; • (ACA T 5tp' + ACE59) 

+\ e ■ [e t ca t Sep' + (e t ce + n T un)56 + n T u56'] + x' d ■ (unse + use')} ds 

+ J l {X' v - (AN T '56) + \ e ■ \NA T 5ip' +[H(exn)|H(wx m)] 56 + M56'} } ds 

(24) 

where we denoted S(a x b) = (a (g> b) — (a ■ b)I, as well as Mv = m x v and 
Nv = n x v, Vv G M 3 ; Moreover, we have 

= D d [£(•)] = D d J(-) + D d G{-) (25) 

where 

(26) 

Finally, we also have 

= D X [L(-)) ■ 5X = J{$K ' An + 6X9 ' ( ETn + nTm ) + 6X 'e • m ) ds ( 2V ) 

For illustration we can consider further that the diameter of a circular cross-section 
is chosen as the design variable, which allows us to express explicitly the result in (23) 
as 



<9n dC, . dC ,. .BA „dA a A, A d\ 

dd = M^- £o) ^M= dmg[E ^ G M ] ' A = — ; (28) 
dm 3D, . dU , ,83 81 8L _ r d 4 rr 

Id = Sd ^ ~ ^ ' ftf = ~dd' E 3 i ' E (^- J=2 ' = l2 (29) 

In order to provide a similar explicit result for directional derivative of the cost func- 
tion, we consider a simple choice given as the beam mass (or equivalently its volume 
for a constant density), 

Ads; A = ^ ( 30 ) 

1 

In such a case the contribution of the cost function to the Kuhn-Tucker optimality 
conditions can be written as 



8 



D d J((j),d)-5d 
D x J((j),d)-5\ 

We can also consider a more complex case of practical interest where the shape op- 
timization is carried out with respect to the beam axis form in the initial configuration. 
A reference configuration is selected in such a case (see Figure 1) and the design vari- 
able is given in terms of the position vector describing the beam initial configuration 
with respect to this reference configuration d = (p (£). The cost function in (25) can 
now be described as 



dA 
~dd 



5dds 



dA 
~dd 



dn 
~2~ 



(31) 
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J(0,d): = Ads= Aj(Z)d£; j(0 



(32) 



6 



In this case all the integrals in (1231) to (1271) must be recomputed with the same 
change of variables like the one presented in ((32l) above typical of the isoparametric 
parent element mapping (e.g. see Zienkiewicz and Taylor OTTO . We also note in 
passing that the derivatives with respect to arc-length coordinate ought to be computed 
by making use of the chain rule 



d ( 
ds 



1 d 



(33) 



For example, the contribution of the cost function to the Kuhn-Tucker optimality 
conditions for such a choice of design variables can be written as 



e 2 



Dd [J{<t>, d)] • 5d = Mjgrd(£) 



gd(Q 
<9£ 



d£ 



(34) 
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3.2 Optimal control 

The optimal control problem studied herein concerns the quasi- static external loading 
sequence which is chosen to bring the structure directly towards an optimal or desired 
final state, which may involve large displacements and rotations. More precisely, we 
study the mechanics problems where introducing the pseudo-time parameter 't' to 
describe a particular loading program is not enough and one also needs to employ the 
control variables v. The latter contributes towards the work of external forces, which 
can be written as 

G ext {u- 54)) := J{5<t> ■ Y Q u}ds (35) 

where F contains the (fixed) distribution of the external loading to be scaled by the 
chosen control. 
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Optimal control can be presented in the following form: find the value of the con- 
trol variables w, such that the final value of the state variables (f> be as close as possible 
to the desired optimal (fixed) values <j> d ; This can be formulated in terms of the con- 
strained minimization of the chosen cost function J(<fi(is), u) which can be written 

as 

J{${v),v)= mm j(<j>(i/*),u*) (36) 

G(cj)(v*);-)=0 4> d 

where the role of the constraint, as indicated by the last expression, is to fix the values 
of the state variables with respect to the chosen value of the control through the given 
set of equilibrium equations. In other words, for a given control we will finally obtain 
the configuration <p the same as the desired state <f> d if the latter verifies equilibrium 
equations, or, in the opposite case simply the solution closest to <f> d which also verifies 
the equilibrium equations. 

In order to remove this constraint and to be able to consider the state variables in- 
dependently from the control variables, we resort to the classical method of Lagrange 
multipliers; Namely, by introducing the Lagrange multipliers A we can rewrite the op- 
timal control problem by making use of the Lagrangian functional L(-) which allows 
to obtain the corresponding form of the unconstrained optimization problem. 



max min L(<f>, v, A) 
vA v(4>,v) 



; L(<j>, u, A) = J{(j>, u) + G(</>, u- A) (37) 

<t>d 



One can readily obtain the Kuhn-Tucker optimality conditions for this problem 
according to 

BL 

= ^-«A:=G(0,MA) (38) 
= « * := A£>. f , + «3**U (39) 

0<p 0<p 0<p 

and 

ou ou ou 
The first of these equations is precisely the weak form of equilibrium equation, whereas 
the second two will provide the basis for computing the control u and the Lagrange 
multipliers A. The explicit form of the equilibrium equation is similar to the one in 
(11) only with the external load term (|35T) . but with the variation of the Lagrange mul- 
tiplier S\ replacing the virtual displacement 8<p. For writing the explicit form of other 
two Kuhn-Tucker equations we choose a particular form of the objective function such 
that 

J(d>, u) = J^H- 0J 2 + ^«IM| 2 } ds (41) 

where cj) d is the desired beam shape and a is a scalar parameter specifying the weighted 
contribution of the chosen control. With such a choice we seek to minimize the "dis- 
tance" between the desired shape and the computed admissible shape (the one which 
satisfies equilibrium) as well as the force or control needed to achieve that state. The 
explicit form of the first term in (|39l) and (|4Q|) can thus be written as 

81 f 

— ■50 = J {5ct> ■ (<t> - <J> d )} ds (42) 
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and 

— ■6v= / {5u-(av)}ds (43) 

Explicit result for the second term of Kuhn-Tucker equations in (l39l) is identical to the 
one in (|24j) potentially modified according to (19) for the case of the follower load. 
Finally, the second term of the Kuhn-Tucker equation in (|40l) can be written as 

QQext 



~iext r 

to Ji 



di 

which concludes with the description of all the problem ingredients. 



(44) 



4 Finite element discrete approximations 

In this section we discuss several important aspects of numerical implementation of 
the presented theory for analysis and design and related issues which arise in numeri- 
cal simulations. 

The analysis part of the problem, i.e. the state variables are represented by using the 
standard isoparametric finite element approximations (e.g. see 11311 ). In particular, this 
implies that the element initial configuration is represented with respect to its parent 
element placed in the natural coordinate space, corresponding to a fixed interval, — 1 = 
6 < f < 6 = +l,by using 



<p (s) = x (0 = }_^N a (0 x a (45) 

a=l 

In (|45T) above x (£) is the position vector field with respect to the reference config- 
uration, x a are nodal values of an element with n en nodes and N a (£ ) are the corre- 
sponding shape functions. The latter can easily be constructed for beams by using the 
Lagrange polynomials, which for an element with n en nodes can be written by using 
the product of monomial expressions 

N a (0= II T^F ^ 

, , , , Stt sfe 
0=1, b=ta 

where £ a , a e [1, n en ] are the nodal values of natural coordinates. 

With isoparametric interpolations one chooses the same shape functions in order 
to approximate the element displacement field, which allows us to construct the finite 
element representation of the element deformed configuration as 

Hen 

cp (0 = $> a (0 ip a (47) 

a=l 

where <p a are the nodal values of the position vector in the deformed configuration. 
The virtual and incremental displacement field are also represented by isoparametric 
finite element interpolations 

S<P (0 = $>« (0 S<p a ; Atp (0 = ]TiV a (0 A<p a (48) 

a=l a=l 
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The latter enables that a new (iterative) guess for the deformed configuration be easily 
obtained with the corresponding additive updates of the nodal values 



Va < ¥>a + Al fia 



The finite element approximation of the incremental displacement field in (1481) 
above, where at each point £ G [£i, £2] tne corresponding value is a linear combination 
of the nodal values, are referred to as the continuum consistent (e.g. see [SdQl), since 
they allow to commute the finite element interpolation and the consistent lineariza- 
tion of nonlinear problem (e.g. lEOlO . Therefore, we also choose the isoparametric 
interpolations for virtual and incremental rotation field with 



56 (£) = J2N a (£) 50 a ; AO (£) = ]TiV a (£) A6 a 



(49) 



a=l 



a=l 



so that the commutativity of the finite element discretization and consistent lineariza- 
tion would also apply to the rotational state variables. The only difference from 
the displacement field concerns the multiplicative updates of the rotation parameters, 
which can be written for any nodal point 'a' as 



A a <- A a exp [AO a ] 



(50) 



In the combined analysis and design procedure proposed herein one must also in- 
terpolate the Lagrange multipliers, which is also done by using the isoparametric in- 
terpolations according to 



A (£) = £> a (0 A " 



Cl=l 



0=1 



(51) 



A e (£) = EJV.(0\ 

0=1 



The corresponding integrals appearing in governing Lagrangian functional in (1211) 
or Kuhn-Tucker optimality conditions in (1231) . (1251) and ((271) are computed by numer- 
ical integration (e.g. Gauss quadrature, see [31]). To illustrate these ideas we further 
state a single element contribution to the analysis part of the governing Lagrangian 
functional in ([22"1) given in the discrete approximation setting by 



G(0 a ,d,A a ) := 



flip 

E 

i=x 



ds 





I 











N a m 



Afe) E(£ ; ) 

n(£,) 1 



T 



m(l) b'(£M-GUA(6)) 



(52) 



where £2 and wi are the abscissas and weights of the chosen numerical integration rule 
(e.g. see [31]) and n ip is the total number of integration points for a single element. 

In order to complete the discretization procedure one must also specify the inter- 
polations of the design variables. If the latter is the thickness or the diameter for the 
chosen case of a circular cross-section, or the element nodal coordinates, it is possible 
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to use again the isoparametric finite element approximations. However, the best results 
are obtained by reducing the number of design variables as opposed to those chosen 
at the element level, by using the concept of a design element. This implies increasing 
the degree of polynomial by employing, for example, Bezier and B-spline curves for 
representation of beam shape and reducing significantly the number of design param- 
eters (e.g. see lfT6l for a detailed discussion of these ideas). What is important to note 
from the standpoint of the simultaneous solution procedure presented further on is that 
the design variable at any point is given as a linear combination of the design element 
interpolation parameters 

d(&) = ]>> a (6) d a (53) 

a=l 

Consequently, the finite or rather design element discretization and consistent lin- 
earization will commute again. This observation was already made earlier for linear 
analysis problem by Chenais and Knopf-Lenoir [3J. With these results on hand we 
can write the discrete approximation of equilibrium equations as 



D\L(<fi, d; A) :=0 
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(54) 



where 'A' denotes the finite element assembly operator. 

By the same token the discrete approximation of the Kuhn-Tucker optimality con- 
ditions in (|23T) can be presented as 



ZVL(0,d;A) :=0 

n e i 
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Similarly, for the simple choice of the objective function in OOl) the discrete ap 
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proximation of the optimality condition in (1251) can be written as 

D d L(<p, d; A) := i-> 
r „ 



dN a (&) . 

ds 




-I T 
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+ B b (Qd b j(Q Wl \ = (56) 



6=1 



In summary, the discrete approximation of the optimal design problem reduces to 
the following set of nonlinear algebraic equations where the unknowns are the nodal 
values of displacements and rotations, the corresponding values of Lagrange mul- 
tipliers and the chosen values of thickness parameters which support Bezier curve 
thickness interpolations, 



r(</>, d; A) :-- 





£int 




r 4> 


:= KA 






dd 


, dV 
~t~ dd 







(57) 



Similar procedure leads to the discrete approximation of the optimal control prob- 
lem. In fact, the first of the nonlinear algebraic equations in (54) is only slightly 
modified to include the control variables in accordance with the expression in (l35l) 



pint 



F i> = 



(58) 



where v are the chosen control parameters. 

For the particular choice of the objective function of the control problem given in 
(|4T|) . we can further obtain the discrete approximation of the optimality condition in 
(l39l) according to 



D <j> L(<f>,v,\) = 



<fi — <fi d + KA = 



(59) 



The discrete approximation of the last optimality condition in (1401) can also be 
written explicitly as 



D v L(<f>, jy, A) = 



au - Fq A 







(60) 



We note in passing that the last two optimality conditions can be combined to elim- 
inate the Lagrange multipliers leading to 



(61) 



The last result can be combined with the equilibrium equations in (1581) providing 
a reduced set of equations with nodal displacements, rotations and control variables 
as the only unknowns. This kind of form is fully equivalent to the arc-length solution 
procedure, which is used for solving the nonlinear mechanics problems in the presence 
of the critical points (e.g. see j[24|. H| or ll23l . among others). One can recognize 
that (1581) is only a particular choice of the supplementary condition which is used to 
stabilize the system. 
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5 Solution procedure 



Two novel solution procedures are developed for solving this class of problems of 
optimal design and optimal control as described next. 

5.1 Diffuse approximation based gradient methods 

The first solution procedure is a sequential one, where one first computes grid values 
of the cost function and then carry out the optimization procedure by employing the 
approximate values interpolated from the grid. It is important to note that all the grid 
values provide the design or control variables along with the corresponding mechani- 
cal state variables of displacements and rotations which must satisfy the weak form of 
equilibrium equation. In other to ensure this requirement, for any grid value of design 
or control variables one also has to solve associated nonlinear problem in structural 
mechanics. 

The main goal of the subsequent procedure is to avoid solving these nonlinear me- 
chanics problems for other but grid values, and simply assume that the interpolated 
values of the cost function will be "sufficiently" admissible with respect to satisfying 
the equilibrium equations. Having relaxed the equilibrium admissibility requirements 
we can pick any convenient approximation of the cost function, which will simplify 
the subsequent computation of the optimal value and thus make it much more efficient. 
These interpolated values of the cost function can be visualized as a surface (yet re- 
ferred to as the response surface) trying to approximate sufficiently well the true cost 
function. The particular method which is used to construct the response surface of 
this kind is the method of diffuse approximations (see E2l or [|2]|). By employing the 
diffuse approximations the approximate value of the cost function J appr is constructed 
as the following quadratic form 

^appr(x) = c + x T b + ^x T Hx (62) 

where c is a constant reference value, b = (£>*) ; hi — 8J Q^ pr is the gradient and H = 

[Hij] ; Hij = d/^T is the Hessian of the approximate cost function. In ((621) above 
variables x should be replaced by either design variables d for the case of an optimal 
design problem or by control variables v in the case when we deal with an optimal 
control problem. 

We further elaborate on this idea for a simple case where only 2 design or control 
variables are used, such that x = (£i,£ 2 ) T - For computational proposes in such a 
case one uses the polynomial approximation typical of diffuse approximation (see 
0) employing the chosen quadratic polynomial basis p(x) and a particular point 
dependent coefficient values a(x) 

J apP r{x) = p T (x)a(x) ; p T (x) = [1, xx, x 2 , xj, x t x 2 , x 2 2 ] ; a(x) = K(x), a 2 (x), . . . , a 6 ( 

(63) 

By comparing the last two expressions one can easily conclude that 
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The approximation of this kind is further fitted to the known, grid values of the cost 
function; J(xj) , % — 1, 2, . . . , n, trying to achieve that the point-dependent coefficients 
remain smooth when passing from one sub-domain to another. This can be stated as 
the following minimization problem: 

1 n 2 

a(x) = arg mm /(a*, x) ; /(a*, x) := - W(x, X;) [J (xj) - p T (x, ; )a*] (65) 

i=i 

where W(x, Xj) are the weighting functions associated with a particular data point 
x, which are constructed by using a window function p(-) based on cubic splines 
according to 

W(x, Xj) = p ( 1 - J ; p(s) : = 1— 3s 2 +2s 3 ; r(x) = max[c?ist(x, x fc )] (66) 

V r ( x ) / fc 

with Xjfc(x), /c = 1, . . . n + 1 (= 7 for the present 2-component case) as the closest 
grid nodes of the given point x. We can see that the weighting functions W(-) in 
(1651) and ((66l) above take a unit value at any of the closest grid nodes x, and vanish 
outside the given domain of influence. While the former assures the continuity of the 
coefficients a(x), the latter ensures that the approximation remains local in character. 
Similar construction can be carried out for higher order problems, which requires an 
increased number of closest neighbors in the list. 

By keeping the chosen point x fixed and considering the coefficients a of diffuse 
approximation as constants, the minimization of /(•) amounts to using the pseudo- 
derivative of diffuse approximation (see li22l ) in order to compute x yielding the min- 
imum of J apP ( x ) according to 

n 

= ^p(x 4 ) ^(x,x l )(J, - p T (x ?; )a) (67) 

i=l 

which allows us to write a set of linear equations 

a(x) = (PWP T )" 1 PWj 
P = [p( Xl ), p(x 2 ) . . . P (x n )] ; j = (J(xi), J(x 2 ), . . . , J(x n )) (68) 
W = diag (W{x 1} x), W(x 2 , x), . . . , W(x n , x)) 

We note in passing that the computed minimum value of J appr does not necessarily 
provide the minimum of the true cost function, which also has to satisfy the equi- 
librium equations; however, for a number of applications this solution can be quite 
acceptable. If the latter is not sufficient, we ought explore an alternative solution pro- 
cedure capable of providing the rigorously admissible value of any computed minima 
of the cost function, by carrying out the simultaneous solution of the cost function 
minimization and the equilibrium equations. The proposed procedure is based on ge- 
netic algorithm as described next. 

5.2 Genetic algorithm based method 

Genetic algorithms belong to the most popular optimization methods nowadays. They 
follow up an analogy of processes that run in nature within the evolution of living 
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organisms over a period of many millions of years. Unlike the classic gradient op- 
timization methods, genetic algorithms operate on so called population which is a 
set of possible solutions, applying the genetic operators, such as cross-over, mutation 
and selection. The principles of genetic algorithm was first proposed by Holland J6J. 
Ever since, the genetic algorithms have reached wide application domain (e.g. see the 
books of Goldberg [Q and Michalewicz [|2TI for extensive review). 

Genetic algorithms in the original form operate on population of so-called chromo- 
somes. These are binary strings which represent possible solutions in a certain way. 
In engineering problems we are usually working with real variables, which in the kind 
of applications described further are optimized values of load control or design vari- 
ables. Adaptation of the genetic algorithm idea to this problem is made possible by 
Storn ll28l by considering the chromosomes as vectors instead of binary strings and 
using differential operators which can affect the distance between the chromosomes. 

In this work we employ an improved version of this kind of algorithm referred to as 
'Simplified Atavistic Differential Evolution', or SADE ([18]). It is shown QUI that 
the algorithm is well suitable for dealing with a fairly large number of variables which 
is of particular interest for the problem on hands. The SADE algorithm is designed 
to explore all possible minima and thus find the global minimum, even for the case 
where the cost function may have very steep gradients and isolated peak values. Only 
a short description of the corresponding procedure is given subsequently (see [7] for a 
more elaborate presentation). 

In the tradition of evolutionary methods, the first step is to generate a starting gen- 
eration of chromosomes by choosing the random values of all state variables along 
with those for control or design. Subsequently we repeat until convergence the cycles 
containing: the creation of a new generation of chromosomes, by mutation, local mu- 
tation or cross-over and the evaluation and selection, which reduces the actual number 
of chromosomes to the initial number. 

In the computations to follow we will work with the population of '10 x n' chromo- 
somes, where 'n' is the total number of unknowns in the problem. This population can 
evolve through the following operations; Mutation: let Xj(g) be the z-th chromosome 
in a generation g, 

*i(g) = (xa(g),x i2 (g),...,Xi n (g)), (69) 

where n is the number of variables of the cost function. If a certain chromosome 
Xj(g) is chosen to be mutated, a random chromosome RP is generated from function 
domain and a new one x fc (g + 1) is computed using the following relation: 

x fe (<? + 1) = x^) + MR(RP - Xi(#)). (70) 

Parameter MR is the constant of algorithm equal to 0.5. Number of new chromosomes 
created by the operator 'mutation' is defined by 'radioactivity', which is another pa- 
rameter of the algorithm, with a constant value set to 0.1 for all our calculations. If 
a certain chromosome is chosen to be locally mutated, all its coordinates have to be 
altered by a random value from a given (usually very small) range. Its aim is to uti- 
lize the near neighborhood of existing chromosomes and make search for improved 
solutions faster. It is useful only for cost functions with steep gradients in the case, 
where near the optimal function value, a small change in value of variable may in- 
troduce a large change in function value. The aim of the cross operator is to create 
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as many new chromosomes as there were in the last generation. The operator creates 
new chromosome Xj (g + 1) according to the following sequential scheme: choose ran- 
domly two chromosomes x q (g) and x r (g), compute their difference vector, multiply it 
by a coefficient CR and add it to the third, also randomly selected chromosome x p (g), 
i.e, 

Xi{g + 1) = x p ( 5 ) + CR(x q (g) - x r ( 5 )). (71) 

Every component exceeding the defined interval is changed to the appropriate 
boundary value of domain. The parameter CR has probably the most important effect 
on the algorithm's behavior; it seems that if the speed up of the convergence is needed, 
this parameter should be set to some lower value (from 0.05 to 0.1). In the opposite 
case, the higher values of this parameter could improve the ability to solve for any 
local minimum. In our computations this value was set to 0.3. 

Selection represents the kernel of the genetics algorithm. The goal is to provide 
a progressive improvement of the whole population, which is achieved by reducing the 
number of "living" chromosomes together with conservation of better ones. Modified 
tournament strategy is used for this purpose: two chromosomes are chosen randomly 
from a population, they are compared and the worse of them is cast off. This conserves 
population diversity thanks to a good chance of survival even for badly performing 
chromosomes. 

It was observed that the SADE algorithm can be quite inefficient in the later stage 
of the analysis, where the solution with all or a large number of components is al- 
most converged. For that reason we have tried to further improve the performance 
by forcing the algorithm to stick better to these converged values. In particular we 
have experimented with a modified form of the cross operator, which contrary to the 
one in (I7TT) will produce the new chromosome by building on top of the best possible 
previous value according to 

Xi(g + 1) = max(x q (g); x r (#)) + SG.CR(x q (g) - x r (#)). (72) 

where SG is the sign change parameter which is supposed to get the correct orientation 
of the increase with respect to the gradient. Moreover, the parameter CR is no more 
constant, but its values are chosen randomly from interval (0, CL), where CL is new 
parameter of the algorithm, with a smaller influence at its behavior than CR. Operator 
local mutation is now switched off and parameter MR in operator mutation is also no 
more constant, but chosen randomly from interval (0, 1). This algorithm modification 
is subsequently referred to as GRADE. 

6 Numerical examples 

In this section we present several illustrative examples dealing with the coupled prob- 
lems of mechanics and either optimal control or optimal design. The computations 
are carried out by using a mechanics model of 2D geometrically exact beam (see 02] 
or H31) developed either within the MATLAB environment for diffuse approximation 
based solution procedure or within the C ++ SADE computer code for genetic algo- 
rithms. 
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6.1 Optimal control of a cantilever structure in the form of letter 
T 
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Figure 2: T letter cantilever: Initial, final and intermediate configurations 



In this example we study the optimal control problem of deploying initially curved 
cantilever beam in the final configuration which takes the form of letter T. See Figure 
|2] for initial and final configurations indicated by thick lines and a number of inter- 
mediate deformed configurations indicated by thin lines. The chosen geometric and 
material properties are as follows: The diameter of the circular curved part and the 
length of the flat part of the cantilever are both equal to 10; the beam cross-section 
is a unit square; the chosen values of Young's and shear moduli are 12000 and 6000, 
respectively. 

The deployment is carried out by applying a vertical force F and a moment M 
at the end of the curved part of the cantilever. In other words the chosen control is 
represented by a vector v = (F, M) T . The desired shape of the cantilever which takes 
the form of letter T corresponds to the values of force F = 40 and moment M = 205. 
The optimal control problem is then defined as follows. The objective or cost function 
is defined by imposing the desired shape only on displacement degrees of freedom 
with 

J(v) = X - Jj4>(v) - <(> d \\ 2 ds (73) 
which is further recast in the discrete approximation setting as 

-, " e ; 2 

J V) = z E E /e ( u » - u ") T ( u » - u ") ^ 

e=0 o=l 

where u*(i/) are computed and are desired nodal displacements. We note in pass- 
ing that no condition is imposed through the cost function on either rotational degrees 
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of freedom or control vector, which nevertheless introduces no difficulties in solving 
this problem. 

The first solution is obtained by the diffuse approximation based gradient method. 
The calculation of the cost function is first carried out for all the 'nodes' of the fol- 
lowing grids: 5 x 5; 10 x 10; 15 x 15 and 20 x 20. The gradient type procedure is 
then started on the grid and, thanks to the smoothness of the diffuse approximation 
base representation of the approximate value of the cost function, converged with no 
difficulty in roughly 20-40 iterations. The iterative values obtained in the gradient 
method computations are for different grids shown in Figured Grid is constructed for 
the following interval of values of force and moment: 

F G (10, 60) ; Me (175,225). 

We note that all different choices of the grid can result in different solutions, since 
none of the solutions of this kind will satisfy the equilibrium equations. 

Quite large difference between known optimal solution and solutions of response 
surface could be explained by different influence of each control variable to value of 
cost function, see Figure HJ All used grids weren't able to describe small influence of 
force F. 

The second solution method used for this problem employs the genetic algorithm 
based computation. In this computation we have used the same admissible intervals 
like in the previous case for the control variables, force and moment. The compu- 
tations are carried out starting from the random values chosen in this interval and 
stopped when the first value of the cost function J < 10~ 7 is found. 

In order to be able to look into statistics, one hundred runs are performed with each 
one converging to the exact solution. Three types of procedures were tried, with either 
tuning on or off the parameters controlling local mutation (LM) and the gradient-type 
acceleration of sign change (SG). The results are presented in Tabled] 



Letter T problem 


Computation 
with SADE 


Number of fitness calls 


Minimum Maximum Mean Value 


LM-on SG-off 
LM-offSG-off 
LM-off SG-on 


240 2860 648.8 
280 1680 625.4 
220 1640 591.2 



Table 1: T-letter cantilever: SADE algorithm performance 

We can see that the best performance is achieved with the simple gradient-like 
modification of SADE genetic algorithm to accelerate convergence in the latest stage. 
The results obtained with this simple gradient-like modification of SADE algorithm 
can further be improved by using described GRADE algorithm, with a special role of 
the radioactivity parameter. See Table [2] 

The second solution procedure applied to the same example is so-called simulta- 
neous solution of mechanics equilibrium and optimal control equations, written ex- 
plicitly in (1581) to (l60l) with value a = 0. The total numb er of unknowns in this case 
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Grid (5 x 5) 

Solution : F = 60.000 
M = 205.26 
38 evaluations of AD 



F approchee 




Grid (10 x 10) 
Solution : F = 59.073 
M = 204.91 
38 evaluations of AD 




Grid (15x15) Grid (20 x 20) 

Solution : F = 51.218 Solution : F = 47.444 
M = 204.95 M = 204.97 

19 evaluations of AD 15 evaluations of AD 

Figure 3: T letter cantilever: Gradient method iterative computation on a grid. 



is equal to 44, with 2 control variables (the force and the moment), 21 components 
of nodal displacements and rotations and 21 Lagrange multiplier. For computational 
convenience, the problem of solving the set of nonlinear algebraic equations in (l58~l) 
to (|60l) is recast as the minimization statement which allows direct application of the 
genetic algorithm with 

= <^ min r T r (75) 




The solution efficiency of the proposed simultaneous procedure depends on the 
chosen upper and lower bounds of the admissible interval and the initial guess for the 
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Whole scale of values. More detailed about value of optimum. 



Figure 4: T letter cantilever: contour of the cost function. 



Letter T problem 


Computation 
with GRADE 


Number of fitness calls 


Minimum Maximum Mean Value 


radioactivity = 0.0 
radioactivity = 0. 1 
radioactivity = 0.2 
radioactivity = 0.3 


180 880 410.0 
240 1200 440.8 
280 1180 512.4 
280 1300 588.4 



Table 2: T letter cantilever: GRADE algorithm performance 

solution. For example, the mechanics state variables are chosen as these featuring in 
the desired beam shape, ip d , and the bounds are controlled by the chosen parameter 
EP according to 

(pe[(l- EP)cp d , (1 + EP)<p d ] (76) 

The equivalent bounds on control variables are obtained from the precious result 
obtained by solving the grid minimization problem which results with the minimum of 
the response surface. Finally, the Lagrange multipliers is solved for from (l59l where 
the adopted values for <p and v are chosen. 

One hundred computations is performed with the indicated bounds on the un- 
knowns. Value of EP parameter is set to 0.0001. Choice of parameters of GRADE 
algorithm were: PR = 30, CL = 2 and 'radioactivity' equal to 0.2. Table [3] summa- 
rizes the statistics of this computation. 





Minimum 


Maximum 


Mean Value 


Standard deviation 


F 


39.957 


40.048 


40.000 


0.0199 


M 


204.86 


205.15 


205.00 


0.0548 


number of evaluations of J(-) 


100740 


298080 


172176 





Table 3: T letter cantilever : solution statistics 



The last part of the analysis carried out in this examples concerns an attempt to 
further increase the efficiency of the simultaneous solution procedure. In that sense, 
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we employ the GRADE version of genetic algorithm with the choice of parameters 
PR = 20, CL = 1 and 'radioactivity' equal to 0.2. A very small value of bounds 
is chosen as well with EP = 0.00001. The computations we have performed with a 
hundred runs of the genetic algorithm can be summarized with the statistics given in 
Table 4. 





Minimal 


Maximal 


Mean Value 


Standard deviation 


F 


39.973 


40.034 


40.000 


0.0135 


M 


204.96 


205.05 


205.00 


0.0192 


number of evaluations of «/(•) 


14720 


201480 


37701 





Table 4: T Letter cantilever : solution statistics 



We can see from Table 4 that the proper choice of the bounds can force the al- 
gorithm to always converge to the same solution. The latter is the consequence of 
using the simultaneous solution procedure which assures that the computed solution 
also satisfies the equilibrium equations. Moreover, the total cost of the simultaneous 
solution procedure can be reduced beyond the one needed for response- surface-based 
approximate solution computations, by either reducing the interval as done herein or 
by making gradient-type modification of this algorithm in order to accelerate the con- 
vergence rate. 

6.2 Optimal control of a cantilever structure in form of letter I 

In the second example we deal with a problem which has a multiple solution and 
its regularized form which should restore the solution uniqueness. To that end, a 
cantilever beam is used very much similar to the one studied in the previous example, 
except for a shorter straight bar with length equal 2. The cantilever is controlled 
by a moment M and a couple of follower forces which follow the rotation of the 
cross-sections to which they are attached. The initial and final configuration, which 
is obtained for a zero couple and a moment M = 205.4 are shown in Figure [51 along 
with a number of intermediate configurations. 

The first computation is performed with the cost function identical to the one in 
(1741) . imposing only the minimum of difference between the desired and the computed 
deformed shape, with no restriction on control variables. The computation is carried 
out by using the GRADE genetic algorithm by starting with random values within the 
selected admissible intervals for the force couple and moment according to 

Fe(0,20); Me (0,230) 

The algorithm performance is illustrated in Tabled 

A number of different solutions have been obtained with different computer runs 
which were performed (see Figure [6]). However, all these solutions remain in fact 
clearly related considering that the applied moment and the force couple play an equiv- 
alent role in controlling the final deformed shape; It can be shown for this particular 
problem that any values of force and moment which satisfy a slightly perturbed ver- 
sion (because of the straight bar flexibility) of the force equilibrium 

F ■ h + M = M = 205.4 
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Figure 5: I letter cantilever: initial, final and intermediate configurations 



Letter I problem 

Number of fitness calls 
Minimum Maximum Mean Value 
180 640 35^8 



Table 5: I letter cantilever: GRADE algorithm performance 

will be admissible solution, thus we have infinitely many solutions for the case where 
only the final shape is controlled by the choice of the cost function. 

In order to eliminate this kind of problem we can perform a regularization of the 
cost function, by requiring not only that the difference between computed and final 
shape be minimized but also that control variables be as small as possible; Namely, 
with a modified form of the cost function 

1 ri el 2 

J h (u a (u), ^) = i EE /e W") - <f (<*.(") - O + ™?*> (77) 

e=0 a=l 

where a is a chosen weighting parameter specifying the contribution of the control. 
We set a very small value a = 10~ 9 and choose the convergence tolerance to 1CT 12 
and carry out the computation for yet a hundred times. Whereas a more stringent 
value of the tolerance requires somewhat larger number of cost function evaluation, 
the result in each of the runs remains always the same, given as 

F = 68.466 ; M = 68.526 

and the found optimal value of the cost function is J = 1.4078631. 10~ 5 . 
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x results of GRADE algorithm 
linear regression for results values: F = 102.63 - 0.49989.M 




moment M 



Figure 6: I letter cantilever: 100 different solutions 



Letter I problem extended 

Number of fitness calls 
Minimum Maximum Mean Value 
820 16320 1758.8 



Table 6: I letter cantilever: GRADE algorithm performance 

6.3 Optimal control of deployment of a multibody system 

The optimal control procedure of the deployment problem of a multibody system is 
studied in this example. In its initial configuration the multibody system consists of 
two flexible component (with EA = 0.1, GA = 0.05 and EI = 10 3 ) each 5 units long 
connected by the revolute joints (see lfi~5l0 to a single stiff component (with EA = 1, 
GA = 0.5 and EI = 10 5 ) of the length equal to 10, which are all placed parallel 
to horizontal axis. In the final deployed configuration the multibody system should 
take the form of a letter B with the stiff component being completely vertical and two 
flexible component considerably bent. The deployment of the system is controlled by 
five control variables, three moments Mi, M 2 and M 3 , a vertical V and a horizontal 
force H. See Figure [7] 

The cost function in this problem is again chosen as the one in|74l which controls 
that the system would find the configuration as close as possible to the desired config- 
uration. The desired configuration of the system corresponds to the values of forces 
H = 0.04, V = -0.05 and moments M 1 = 0.782, M 2 = -0.792 and M 3 = 0.792. 
The solution is computed by using the SADE and GRADE genetics algorithms and 
starting with the random choice in the interval of interest defined as 

H E (0.025; 0.05), V E (-0.06; -0.035), M x E (0.6; 0.9), 
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EA 1= 1.0 EA 2 =0.1 
cGAi=1.0 cGA 2 =5.0 
Eli =1.05 EI 2 = 1000.0 



Final shape 




Figure 7: Multibody system deployment: initial, final and intermediate configurations. 

M 2 G (-0.9; -0.65), M 3 G (0.6; 0.85). 

The solution of the problem is typically more difficult to obtain with an increase in 
the number of control variables, one of the reasons being a more irregular form of the 
cost function. In that sense,we refer to illustrative representation of the cost function 
contours in different subspaces of control variables as shown in Figures [8] and |9] 

The convergence tolerance on cost function is chosen equal to 10~ 6 . The SADE 
algorithm performance for the simplest choice of algorithm parameters is presented 
in Tableland the GRADE algorithm performance for modified value of radioactivity 
parameter is presented in Tabled 



Letter B problem 


Computation 
with SADE 


Number of fitness calls 


Minimum Maximum Mean Value 


LM-on SG-off 
LM-offSG-off 
LM-off SG-on 


2600 165800 46887.5 
2350 177150 42494.0 
2400 177850 34612.1 



Table 7: Results of SADE algorithm for 5-variable control problem 

One can notice the order of magnitude of increase in cost function evaluation, 
which is brought about by a more elaborate form of the cost function (see Figures 
[8] and [9]). However, the latter is not the only reason. In this particular problem the role 
of moments in the list of control variables is much more important than the role of the 
horizontal and vertical forces in bringing the system in the desired shape. This affects 
the conditioning of the equations to be solved and the slow convergence rate of the 
complete system is in reality only the slow convergence of a single or a couple of con- 
trol components. The latter is illustrated in Figure \\0[ where we provide the graphic 
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V — Mi subspace V — Mi subspace 

Figure 8: Multibody system deployment: contours of the cost function in different 
subspaces. 

representation of iterative values for computed chromosomes, where every chromo- 
some is represented by a continuous line. We can note that the population of optimal 
values of moments converges much more quickly than the force values which seek 
a large number of iteration in order to stabilize. Another point worthy of further ex- 
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Mi — M 2 subspace 



Mi — M 2 subspace 




M 2 — M 3 subspace M 2 — M 3 subspace 

Figure 9: Multibody system deployment: contours of the cost function in different 
subspaces. 



ploration is the best way to accelerate the convergence rate in the final computational 
phase. 
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Letter B problem 


Computation 
with GRADE 


Number of fitness calls 


Minimum Maximum Mean Value 


radioactivity = 0.0 
radioactivity = 0.1 
radioactivity = 0.2 
radioactivity = 0.3 


1500 114050 23862.5 
1900 117850 20632.0 
3050 122550 34520.0 



Table 8: Results of GRADE algorithm for 5D task 




G1:J=0.0091823 G4:J=0.0071823 




G26:J=1.3302e-05 G43:J=1.1839e-05 

Figure 10: Multibody system deployment: convergence of iterative chromosome pop- 
ulations 



6.4 Optimal design of shear deformable cantilever 

In this last example we study an optimal design problem which considers the thickness 
optimization of a shear deformable cantilever beam, shown in Figure [TT1 

The beam axis in the initial configuration of the cantilever and the thickness is 
considered as the variable to be chosen in order to assure the optimal design speci- 
fied by a cost function. In the setting of discrete approximation, we choose 4 beam 
elements each with a constant thickness hi, which results with 4 design variables 
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h 1 


h 2 


h 3 


h 4 

S 


' - 










- 


E = 75000 










G = 50000 










c = 5/6 










p = 1/30 


b = 30 








imposed mass: 


Mo = 30000 









-100 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

-200 200 400 600 800 1000 1200 

Figure 1 1 : Shear deformable cantilever beam optimal design : initial and deformed 
shapes 



d = h = (hi, h 2 , h 3 , h^j. The beam mechanical and geometric properties are: 
Young's modulus E = 75000, shear modulus G = 50000, rectangular cross section 
b x hi with width b = 30 and mass density p = 1/30. The latter is needed for comput- 
ing the total mass of the beam M = J L pbh(s)ds to be used as the corresponding lim- 
itation on the computed solution assuring reasonable values of the optimal thickness 
under the free-end vertical force F = 1000. In order to assure a meaningful result the 
computations are performed under chosen value of mass limitation is M = 30000. 
Other limitations are also placed on the admissible values of the thickness for each 
element. 

The first computation is performed by using the diffuse approximation based re- 
sponse function and the sequential solution procedure. The cost function is selected 
as the shear energy of the beam and problem is cast as maximization of the shear 
energy, with 

J(<p(d)) = max J(£*(d*)) ; J(<p m ) = [ iGA^ds (78) 

where 7 is the shear strain component. The bounds on thickness values are chosen 
as shown in Tabled The diffuse approximation computations on the grid are started 



Thickness 






h 


hi 


Min 


30 


30 


15 


15 


Max 


60 


60 


35 


35 



Table 9: Shear deformable cantilever optimal design : thickness admissible values 

from an initial guesses for thickness h = (55,50,30,20). It took 11 iterations to 
converge to solution given as 

h = (45.094, 40.074, 19.832, 15.000) 
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The corresponding value of shear energy for this solution is J appr = 16.3182; we recall 
it is only an approximate solution, since the computed value does not correspond to 
any of the grid nodes. 

The same solution is next sought by using GRADE algorithm; The chosen values 
for GRADE parameters are: PR = 10, CL = 1.0 and 'radioactivity' = 0.2. The ge- 
netic algorithm is executed 100 times, leading to the computational statistics reported 
in Table C0l 





Minimum 


Maximum 


Mean Value 


nb. of comput. J(-) 


120 


3400 


674.8 



Table 10: Shear deformable cantilever optimal design : computation statistics 



The algorithm yielded two different solutions, both essentially imposed by the cho- 
sen bounds; Namely, out of 100 runs, 57 converged to h = (60, 30, 15, 15), with the 
corresponding value of J = 17.974856, whereas 43 settled for h = (30, 60, 15, 15) 
with a very close value of J = 17.926995. Hence, each of two solutions leads to an 
improved value of the cost function. 

The second part of this example is a slightly modified version of the first one, in 
the sense that the mechanics part of the problem is kept the same and only a new 
cost function is defined seeking to minimize the Euclidean norm of the computed 
displacement vector, i.e. 

J(<p(d)) = max J(£*(d*)) ; J(cp*) = ~ / ||p - x|| 2 d S (79) 
G(v?*(d*),-)=o 2 J L 

Such a choice of the cost function is made for being well known to result with a 
well-conditioned system expressing optimality conditions. Indeed, the same type of 
sequential solution procedure using diffuse approximation of cost function now needs 
only a few iterations to find the converged solution, starting from a number initial 
guesses. The final solution value is given as 

h = (42.579, 35.480, 26.941, 15.000) 



In the final stage of this computation we recompute the solution of this problem by 
using the genetic algorithm. The admissible value of the last element thickness is also 
slightly modified by reducing the lower bound to 5 (instead of 15) and higher bound 
to 25 (instead of 35) in order to avoid the optimal value which is restricted by a bound. 
The first solution to this problem is obtained by using again the sequential procedure, 
where the GRADE genetic algorithm is employed at the last stage. The computed 
value of the displacement vector norm for found solution is 623808 and mass M is 
30062. The computations are carried out a hundred times starting from random initial 
values. The statistics of these computations are given in Table [TT1 

The same kind of problem is now repeated by using the simultaneous solution pro- 
cedure, where all the optimality condition are treated as equal and solved simultane- 
ously resulting with 4 thickness variables, 15 displacement and rotation components 
and as many Lagrange multipliers as unknowns. The latter, in fact, is eliminated prior 
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Minimum 


Maximum 


Mean Value 


Standard deviation 


7 

hi 


43.772 


43.807 


43.790 


0.0094 


7 

h 2 


35.914 


35.949 


35.932 


0.0088 


h 

ri3 


Of 'J 1 'J 


Of 'ZAf 
ZD. .540 


zo. Jzo 


U.UUoz 




14 1 84 


14 210 


14 197 

± ■~r . J. / 


0064 


nb. of comput. «/(•) 


1440 


9960 


3497 





Table 1 1 : Shear deformable cantilever optimal design : computation statistics 

to solution by making use of optimality condition in (f57|) ?. The chosen upper and 
lower bounds of the admissible interval are chosen as 

Lpe[(l-EP)^,(l+EP)<f?} (80) 

where the guess for the displacement ip v is obtained by solving mechanics problem 
with the values of thickness parameters given in Table ITTl The limitation on total mass 
is added to the cost function. The choice of GRADE algorithm parameters is given 
as PR = 20, CL = 2 and 'radioactivity' equal to 0.1. The computation is stopped 
with a fairly loose tolerance, which allows to accelerate the algorithm convergence but 
does not always lead a unique solution. Yet, the results in Table [Tzl show that standard 
deviation indeed remains small, or that the solution is practically unique. 





Minimum 


Maximum 


Mean Value 


Standard deviation 


hi 


43.782 


43.794 


43.789 


0.0026 


h 2 


35.925 


35.935 


35.930 


0.0021 


h 


26.315 


26.324 


26.319 


0.0019 


h^ 


14.197 


14.202 


14.200 


0.0010 


nb. of comput. r T r 


111340 


968240 


313006 





Table 12: Shear deformable cantilever optimal design : simultaneous computation 
statistics 
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8 Conclusions 

The approach advocated herein for dealing with a coupled problem of nonlinear struc- 
tural mechanics and optimal design or optimal control, which implies bringing all the 
optimality conditions at the same level and treating all variables as independent rather 
than considering equilibrium equations as a mere constraint and state variables as de- 
pendent on design or control variables, is fairly unorthodox and rather unexplored. 
For a number of applications the proposed approach can have a great potential. In 
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particular, the problems of interest to this work concern the large displacements and 
rotations of a structural systems. The key ingredient of such an approach pertains to 
geometrically exact formulation of a nonlinear structural mechanics problem, which 
makes dealing with nonlinearity description or devising the solution schemes much 
easier then for any other model of this kind. The model problem of the geometrically 
exact beam explored in detail herein is not the only one available in this classl; We 
refer to work in [9 J for shells or to ifTOll for 3D solids, with both models sharing the 
same configuration space for mechanics variables as 3D beam. The latter also allows 
to directly exploit the presented formulation and the solution procedures of a coupled 
problem of nonlinear mechanics, for either shells or 3D solids, and optimal control or 
optimal design. 

Two different solution procedures are presented herein; the first one, which exploits 
the response surface representation of the true cost function followed by a gradient 
type solution step, leads to only an approximate solution. Although the quality of 
such a solution can always be improved by further refining the grid which serves 
to construct the response surface, the exact solution is never computed unless the 
minimum corresponds to one of the grid points. The second solution procedure, which 
solves simultaneously the optimality conditions and nonlinear mechanics equilibrium 
equations, does deliver the exact solution, although often only after the appropriate 
care is taken to choose the sufficiently close initial guess and to select the admissible 
intervals of all variables accordingly. Probably the best method in that sense is the 
combination of sequential and simultaneous procedure, where the first serves to reduce 
as much as possible the admissible interval and provide the best initial guess, whereas 
the second furnishes the exact solution. 

A number of further improvements can be made for the proposed methods of this 
kind, to help increase both the convergence rates and accuracy of the computed solu- 
tion. One has to remember that even only mechanics component equilibrium equations 
for a problem of this kind can be very difficult to solve, since the large difference in 
stiffness in different modes (e.g. bending versus stretching) can result with a poorly 
conditioned set of equations. Adding the optimality conditions on top only further 
increases the difficulty. Finding the best possible way to deal with this problem is 
certainly worthy of further explorations. 
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